gisR

# Lesen der Gebäudedaten
gwr <- read.csv(paste0("https://www.web.statistik.zh.ch/ogd/daten/",
                       "ressourcen/KTZH_00002022_00004064.csv")
                ) |>
  dplyr::rename_all(tolower)


# Umwandeln in ein sf-Objekt -> Geodatenobjekt
gwr_sf <- sf::st_as_sf(
  gwr, 
  coords = c("e.gebaeudekoordinate", "n.gebaeudekoordinate"), 
  crs = 2056
)

# interaktives Darstellen
mapview::mapview(head(gwr_sf, 100))

gemeinden <- sf::st_read(
  paste0("https://maps.zh.ch/wfs/OGDZHWFS?SERVICE=WFS&VERSION=2.0.0&",
         "Request=getfeature&TYPENAME=ms:ogd-0095_arv_basis_up_gemeinden_seen_f&",
         "outputformat=geojson")
  ) |>
  dplyr::rename_all(tolower)

# Matchen der Gemeindeinformation an die GWR-Daten
# Die Geometrie wird hier entfernt, da die Datenanalyse danach 
# viel schneller läuft
gwr_sf_gem <- gwr_sf |>
  dplyr::select(
    gebaeudekategorie_code, 
    gebaeudekategorie_bezeichnung
    ) |>
  sf::st_join(gemeinden) |>
  sf::st_drop_geometry()


# Berechnung der Anteile der Gebäude ohne Wohnnutzung 
# am gesamten Gebäudebestand einer Gemeinde
geb_ohne_wohnnutzung_pro_gem <- gwr_sf_gem |>
  dplyr::group_by(
    gebaeudekategorie_code, 
    gebaeudekategorie_bezeichnung, 
    bfs, 
    gemeindename
  ) |>
  dplyr::summarise(anzahl = dplyr::n()) |>
  dplyr::ungroup() |>
  dplyr::group_by(bfs, gemeindename) |>
  dplyr::mutate(anteil = round(anzahl/sum(anzahl)*100, 2)) |>
  dplyr::ungroup() |>
  dplyr::filter(gebaeudekategorie_code == 1060)

# Hinzufügen der Gemeindepolygone
# Herausfiltern der Seen und fixen der Polygone
geb_ohne_wohnutzung_sf <- gemeinden  |>
  dplyr::select(bfs) |>
  dplyr::left_join(geb_ohne_wohnnutzung_pro_gem, by = "bfs") |>
  dplyr::filter(bfs != 0) |>
  sf::st_make_valid() |>
  # tmap vewendet die erste Spalte als Bezeichner, 
  #deshalb hier die Änderung in der Spaltenanordnung
  dplyr::select(gemeindename, dplyr::everything())


geb_ohne_wohnnutzung_sf_generalized <- geb_ohne_wohnutzung_sf |>
  rmapshaper::ms_simplify(keep = 0.005)

## Nur für Präsentation!
geb_ohne_wohnnutzung_sf_generalized <- readRDS(here::here("extdata/data.RDS"))

# Statisches Darstellung
tmap::tm_shape(geb_ohne_wohnnutzung_sf_generalized) +
  tmap::tm_polygons("anteil") +
  tmap::tm_compass() +
  tmap::tm_scale_bar()

## Nur für Präsentation!
geb_ohne_wohnnutzung_sf_generalized <- readRDS(here::here("extdata/data.RDS"))

# Setzt tmap auf interaktiv
tmap::tmap_mode("view")

# Interaktive Darstellung
tmap::tm_shape(geb_ohne_wohnnutzung_sf_generalized) +
    tmap::tm_polygons("anteil")

stac_source <- rstac::stac(
  "https://planetarycomputer.microsoft.com/api/stac/v1",
  force_version = "1.0.0"
)
collections_query <- stac_source |>
  rstac::collections()

rstac::get_request(collections_query)

# BBOX definieren
ktzh_bbox <- sf::read_sf("geodata/Gemeindegrenzen_-OGD.gpkg") |> 
  sf::st_transform("EPSG:4326") |> 
  sf::st_make_valid() |> 
  sf::st_union() |> 
  sf::st_bbox() |> 
  format_bbox()
  

# STAC Query bauen
stac_query <- rstac::stac_search(
  q = stac_source,
  collections = "sentinel-2-l2a",
  datetime = "2023-06-10/2023-06-14",
  bbox = ktzh_bbox
) |>
  rstac::ext_filter(
    `eo:cloud_cover` <= 20)

wfs_zh <- "https://maps.zh.ch/wfs/OGDZHWFS"

url <- httr::parse_url(wfs_zh)
url$query <- list(service = "WFS",
                  version = "2.0.0",
                  request = "GetFeature",
                  typename = "ms:ogd-0095_arv_basis_up_bezirke_f",
                  outputformat = "geojson")

request <- httr::build_url(url)
bezirke <- sf::read_sf(request)

# Save Bezirke Layer
sf::write_sf(bezirke, "geodata/Bezirke.gpkg")

# calculate NDVI using the red (band 1) and nir (band 4) bands
sent2_ndvi <- (RGBI_zh[[4]] - RGBI_zh[[1]]) / (RGBI_zh[[4]] + RGBI_zh[[1]])

sent2_ndvi |> terra::writeRaster("geodata/ndvi_terra.tif")

# Load Polygon Layer
LWNutz_Diels <- sf::read_sf("geodata/LWNutz_Dielsdorf.gpkg") |>
  dplyr::mutate(area = sf::st_area(geom))

# Set aller Flächen, die erst am 15.06. sollten geschnitten werden:
BF_Fl <- LWNutz_Diels |>
  dplyr::filter(harvest_date == "15.06.") |>
  # Drop die kleinsten Flächen unter 500 m^2
  dplyr::filter(area > units::set_units(500, "m^2"))

# Setzt tmap auf interaktiv
border <- c(2669379, 1252226, 2684597, 1269716)
breaks <- c(-0.4, -0.2, 0, 0.2, 0.4, 1)

tmap::tm_shape(BF_Fl, bbox = border) +
tmap::tm_polygons("mean_ndvi", style = "fixed", breaks = breaks) +
tmap::tm_shape(BF_Fl |> dplyr::filter(mean_ndvi <= 0.4), bbox = border, name = "NDVI <= 0.4 ") +
tmap::tm_polygons("mean_ndvi", style = "fixed", breaks = breaks, legend.show = FALSE)

BF_Fl |> dplyr::filter(`mean_ndvi` <= 0.4) |> 
  dplyr::select(blw_name, region, mean_ndvi, area) |>
  sf::st_drop_geometry() |> head(10) |> 
  knitr::kable() |> kableExtra::kable_styling(font_size = 22)